REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 

regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 


2.  REPORT  TYPE 

New  Reprint 


4.  TITLE  AND  SUBTITLE 

Crystal  structure  prediction  for  cyclotrimethylene  trinitramine 
(RDX)  from  ?rst  principles 


3.  DATES  COVERED  (From  -  To) 


5a.  CONTRACT  NUMBER 

W91  INF-09-1-0397 


5b.  GRANT  NUMBER 


6.  AUTHORS 

Ratal  Podeszwa,  Betsy  M.  Rice,  and  Krzysztof  Szalewicz 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

University  of  Delaware 
Vice  Provost  for  Research 
University  of  Delaware 

Newark,  DE  19716  - 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND 
ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

53252-CH.2 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  public  release;  federal  purpose  rights 


13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


14.  ABSTRACT 

Crystal  structure  prediction  and  molecular  dynamics  methods  were  applied  to  the 

cyclotrimethylene  trinitramine  (RDX)  crystal  to  explore  the  stability  rankings  of  various  polymorphs  using  a 
recently  developed  nonempirical  potential  energy  function  that  describes  the  RDX  dimer  interactions.  The  energies 
of  500  high-density  structures  resulting  from  molecular  packing  were  minimized  and  the  14  lowest-energy 
structures  were  subjected  to  isothermal-isostress  molecular  dynamics  (NsT-MD)  simulations.  For  both  crystal 


15.  SUBJECT  TERMS 

symmetry-adapted  perturbation  theory,  crystal  structure 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

15.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF  PAGES 

Krzysztof  Szalewicz 

UU 

UU 

UU 

UU 

19b.  TELEPHONE  NUMBER 

302-831-6579 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39. 1 8 


Report  Title 

Crystal  structure  prediction  for  cyclotrimethylene  trinitramine  (RDX)  from  ?rst  principles 

ABSTRACT 

Crystal  structure  prediction  and  molecular  dynamics  methods  were  applied  to  the 

cyclotrimethylene  trinitramine  (RDX)  crystal  to  explore  the  stability  rankings  of  various  polymorphs  using  a  recently 
developed  nonempirical  potential  energy  function  that  describes  the  RDX  dimer  interactions.  The  energies  of  500 
high-density  structures  resulting  from  molecular  packing  were  minimized  and  the  14  lowest-energy  structures  were 
subjected  to  isothermal— isostress  molecular  dynamics  (NsT-MD)  simulations.  For  both  crystal  structure  prediction 
methods  and  molecular  dynamics  simulations,  the  lowest-energy  polymorph  corresponded  to  the  experimental 
structure;  furthermore,  the  lattice  energy  of  this  polymorph  was  lower  than  that  of  the  other  polymorphs  by  at  least  1 . 1 
kcal/mol. 

Crystal  parameters  and  densities  of  the  low-energy  crystal  produced  by  the  NsT-MD  simulations  matched  those  of  the 
experimental  crystal  to  within  1%  of  density  and  cell  edge  lengths  and  0.01 1  of  the  cell  angle.  The  arrangement  of  the 
molecules  within  the  time-averaged  unit  cell  were  in  equally  outstanding 

agreement  with  experiment,  with  the  largest  deviation  of  the  location  of  the  molecular  mass  centers  being  less  than  0.07 
A  and  the  largest  deviation  in  molecular  orientation  being  less  than  2.81.  NsT-MD  simulations  were  also  used  to 
calculate  crystallographic  parameters  as  functions  of  temperature  and  pressure  and  the  results  were  in  a  reasonable 
agreement  with  experiment. 


REPORT  DOCUMENTATION  PAGE  (SF298) 
(Continuation  Sheet) 


Continuation  for  Block  13 


ARO  Report  Number  53252. 2-CH 

Crystal  structure  prediction  for  cyclotrimethylene 


Block  13:  Supplementary  Note 

©2009  RCS.  Published  in  Physical  Chem  Chem  Physics,  Vol.  11,5512,  (2009),  (5512).  DoD  Components  reserve  a 
royalty-free,  nonexclusive  and  irrevocable  right  to  reproduce,  publish,  or  otherwise  use  the  work  for  Federal  purposes,  and  to 
authroize  others  to  do  so  (DODGARS  §32.36).  The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the 
author(s)  and  should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by 
other  documentation. 

Approved  for  public  release;  federal  purpose  rights 


PAPER 


www.rsc.org/pccp  |  Physical  Chemistry  Chemical  Physics 


Crystal  structure  prediction  for  cyclotrimethylene  trinitramine  (RDX) 
from  first  principlesf 

Rafal  Podeszwa,*"6  Betsy  M.  Rice'  and  Krzysztof  Szalewicz* 

Received  4th  February  2009,  Accepted  25th  March  2009 
First  published  as  an  Advance  Article  on  the  web  28th  April  2009 
DOI:  10.1039/b902015b 


Crystal  structure  prediction  and  molecular  dynamics  methods  were  applied  to  the 
cyclotrimethylene  trinitramine  (RDX)  crystal  to  explore  the  stability  rankings  of  various 
polymorphs  using  a  recently  developed  nonempirical  potential  energy  function  that  describes  the 
RDX  dimer  interactions.  The  energies  of  500  high-density  structures  resulting  front  molecular 
packing  were  minimized  and  the  14  lowest-energy  structures  were  subjected  to 
isothermal-isostress  molecular  dynamics  (NsT-MD)  simulations.  For  both  crystal  structure 
prediction  methods  and  molecular  dynamics  simulations,  the  lowest-energy  polymorph 
corresponded  to  the  experimental  structure;  furthermore,  the  lattice  energy  of  this  polymorph  was 
lower  than  that  of  the  other  polymorphs  by  at  least  1.1  kcal  mol-1.  Crystal  parameters  and 
densities  of  the  low-energy  crystal  produced  by  the  NsT-MD  simulations  matched  those  of  the 
experimental  crystal  to  within  1%  of  density  and  cell  edge  lengths  and  0.01°  of  the  cell  angle.  The 
arrangement  of  the  molecules  within  the  time-averaged  unit  cell  were  in  equally  outstanding 
agreement  with  experiment,  with  the  largest  deviation  of  the  location  of  the  molecular  mass 
centers  being  less  than  0.07  A  and  the  largest  deviation  in  molecular  orientation  being  less  than  2.8°. 
NsT-MD  simulations  were  also  used  to  calculate  crystallographic  parameters  as  functions  of 
temperature  and  pressure  and  the  results  were  in  a  reasonable  agreement  with  experiment. 


1.  Introduction 

Prediction  of  a  crystal  structure  using  only  information  about 
the  atomic  arrangement  of  a  constituent  molecule  was  once 
considered  impossible.1-4  While  the  structure  of  various  con¬ 
formations  of  individual  molecules  can  be  readily  determined 
on  current  computational  architectures  using  highly-correlated 
quantum  mechanical  methods,  there  are  substantial  challenges 
that  preclude  the  prediction  of  the  crystalline  state  structures 
of  even  fairly  small  organic  molecules.  The  issue  lies  with 
properly  representing  the  binding  forces  within  molecular 
crystals:  van  der  Waals  interactions.  Such  forces  can  be 
modelled  ab  initio  by  well-established  methods  of  quantum 
chemistry.  However,  the  high  computational  cost  of  such 
approaches  has  prevented  extensive  applications  of  these 
methods  to  molecular  crystals.  The  only  computationally 
feasible  quantum  mechanics-based  approach  for  treating 
molecular  crystals  is  density  functional  theory  (DFT). 
Unfortunately,  DFT  has  been  shown  to  be  inaccurate  for 
systems  in  which  the  dispersion  components  of  van  der  Waals 
forces  are  dominant  (see,  e.g.,  ref.  5).  This  problem  extends  to 
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DFT  calculations  for  molecular  crystals  with  periodic  bound¬ 
ary  conditions.6’7  Empirical  models  using  simple  atom-atom 
potentials  can  be  used  to  predict  crystal  structures,  but  rely  on 
limited  experimental  information  (possibly  supplemented  by 
some  quantum  mechanical  calculations)  for  parameterization. 
Consequently,  atomistic  simulations  using  such  models  often 
produce  unsatisfactory  results  in  describing  systems  or  condi¬ 
tions  outside  of  those  included  in  the  fitting  data.  All  of  these 
factors  have  made  the  prediction  of  crystal  structure 
using  information  from  only  the  molecular  structure  rather 
unsuccessful.1-4,8 

Recently,  we  demonstrated9,10  that  potential  energy  func¬ 
tions  fitted  to  calculations  using  a  novel  method,  symmetry- 
adapted  perturbation  theory  based  on  the  density-functional 
description  of  the  monomers  [SAPT(DFT)], 11-13  provide  suffi¬ 
cient  accuracy  and  numerical  efficiency  for  predicting  the 
crystal  structure  properties  of  the  molecular  crystal  of  cyclo¬ 
trimethylene  trinitramine  (RDX).  Also,  for  the  benzene  crystal 
at  the  experimental  geometry,  the  SAPT(DFT)  interaction 
energies,14  together  with  three-body  SAPT(DFT)  non-additive 
corrections,15  yielded  the  lattice  energy  accurate  to  within  a 
few  percent  of  the  experimental  values,9  significantly  better 
than  other  methods  based  on  first  principles.16,17  The  work  of 
ref.  9  and  10  was  based  on  a  six-dimensional  potential  fitted  to 
over  1100  single-point  calculations  for  the  RDX  dimer  using 
the  SAPT(DFT)  method.  We  first  used  this  potential  in 
isothermal-isostress  molecular  dynamics  (NsT-MD)  simula¬ 
tions  of  crystalline  RDX  at  ambient  conditions  for  the 
experimentally  observed  structure.10  We  showed  that  the 
crystallographic  parameters  and  molecular  configurations 
were  in  excellent  agreement  with  measured  values.  In  the 
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current  study,  which  extends  and  supplements  the  work 
described  in  ref.  9,  we  use  the  SAPT(DFT)  potential  of  the 
RDX  dimer  in  NsT-MD  simulations  of  several  hypothetical 
RDX  crystal  structures,  where  the  initial  configurations  are 
obtained  using  crystal  structure  prediction  methods.  This 
work  will  explore  several  structural  phases  of  RDX  other  than 
those  observed  experimentally,  and  the  role  of  thermal  effects 
on  the  stability  rankings.  Additionally,  the  temperature  and 
pressure  dependencies  of  crystallographic  parameters  for  the 
experimentally  observed  polymorph  will  be  determined  from 
NsT-MD  simulations. 

Simultaneously  with  our  work,910  the  goal  of  predicting 
crystal  structures  was  achieved  by  two  other  groups. 
The  method  of  Misquitta  et  a/.,18  successfully  applied  in  a 
prediction  of  the  lowest-energy  structure  of  the  crystal  of 
C6Br2ClFH2,  also  uses  SAPT(DFT)  and  can  be  viewed 
as  a  simplification  of  the  approach  used  by  us.  The  main 
differences  are  in  the  treatment  of  the  induction  and  dispersion 
energies  and  of  their  exchange  counterparts.  Whereas  we 
compute  all  these  terms  for  the  whole  set  of  dimer  configura¬ 
tions,  Misquitta  et  al.  approximate  them  by  a  damped  asymp¬ 
totic  expansion.  There  are  also  differences  in  the  generation 
and  ranking  of  polymorphs.  In  both  cases,  a  large  number  of 
polymorph  candidates  is  first  generated,  then  a  subset 
of  polymorphs  is  subjected  to  lattice  energy  minimizations 
(a  0  K  simulation).  In  our  study,  the  minimizations  are 
followed  by  NsT-MD  simulations  at  ambient  conditions,  thus 
introducing  thermal  effects.  We  will  further  elaborate  on  these 
calculations  below. 

A  very  different  approach  was  developed  by  Neumann 
et  al.19  These  authors  first  fitted  a  simple  atom-atom  potential 
function  for  a  given  system  to  lattice  energies  computed  using 
the  DFT  +  D  (DFT  plus  asymptotic  dispersion)  method 
of  Neumann  and  Perrin20  and  applied  it  in  lattice  energy 
minimization.  For  a  number  of  lowest-energy  polymorphs 
predicted  by  this  potential,  they  perform  further  optimi¬ 
zations  with  “on-the-fly”  DFT  +  D  calculations.  This  approach 
works  reasonably  well,  despite  the  known  problems  of  the 
standard  DFT  method  in  computing  interaction  energies 
which  cannot  be  solved  just  by  a  simple  addition  of  true 
dispersion  energies.5  To  overcome  this  problem,  Neumann 
and  Perrin20  introduced  a  “damping"  function  in  the  disper¬ 
sion  component  with  parameters  optimized  to  reproduce 
experimental  lattice  parameters  for  about  30  crystals  of  small 
molecules.  Thus,  this  DFT  +  D  method  includes  a  significant 
degree  of  empiricism,  in  contrast  to  our  SAPT(DFT)-based 
approach  which  is  completely  nonempirical. 

Computational  details 

Crystal  structure  predictions  for  RDX 

Our  first  goal  is  to  generate  probable  RDX  crystal  structures 
using  information  only  about  the  molecular  structure.  Several 
approaches  have  been  developed  to  generate  and  rank  probable 
crystal  structures.21  We  used  a  method  developed  specifically 
for  molecular  energetic  materials22,23  which  should  be  appro¬ 
priate  for  the  system  under  study  here.  In  the  first  step  of 
this  procedure,  a  three-dimensional  molecular  model  of  the 


constituent  molecule  (monomer)  is  generated.  Here,  we  used 
the  rigid  monomer’s  structure  corresponding  to  that  measured 
in  the  crystal  at  ambient  conditions.24  Next,  crystals  of 
different  symmetries  are  created  using  this  monomer’s  struc¬ 
ture  and  the  MOLPAK  (molecular  packing)  software  suite.22 
MOLPAK  produces  crystals  corresponding  to  molecular 
coordination  geometries  observed  in  the  most  common  space 
groups  of  molecular  crystals;  descriptions  of  coordina¬ 
tion  geometries  are  given  in  ref.  22  and  25.  The  version  of 
MOLPAK  used  in  this  study25  samples  51  coordination 
geometries  corresponding  to  the  triclinic  (PI,  Pi),  monoclinic 
(P2U  Ph/c,  Cc ,  C2,  C2/c,  Pc,  Pile,  P2i/m,  P2/m,  P2,  Pm, 
P2/m),  and  orthorhombic  (Z  =  4,  P2{2\2,  P212121,  Pca2\, 
Pna2 1 ,  Pnn2,  Pba2,  Pnc2,  P222\,  Pmn2\,  Pma2,  Z  =  8,  Pbcn, 
Pbca)  space  groups.  These  coordination  geometries  were 
found  via  a  crystal  database  search  as  described  in  ref.  22 
and  Appendix  1  of  ref.  25.  For  each  hypothetical  crystal 
generated,  an  initial  arrangement  of  the  atoms  within  the  unit 
cell  is  obtained  by  first  orienting  the  monomer  in  a  Cartesian 
coordinate  system  whose  origin  is  located  at  the  centroid  of  the 
monomer.  The  remaining  symmetry-equivalent  entities  are 
then  generated  using  the  appropriate  space  group  symmetry 
operations  and  coordinate  geometry  definition  as  described  in 
ref.  22.  Sampling  of  possible  orientations  of  the  molecules 
within  each  coordination  geometry  is  performed  by  a  systematic 
angular  sweep  of  the  range  ±  90°  in  10°  steps.  This  sampling 
results  in  the  generation  of  6,859  (193)  hypothetical  crystal 
structures  per  coordination  geometry.  The  intermonomer 
separations  were  then  varied  for  each  hypothetical  structure 
until  a  minimum-volume  packing  based  on  a  repulsion 
criterion  is  reached  as  described  in  ref.  22  and  Appendix  1  of 
ref.  25.  In  this  way,  the  MOLPAK  step  generated  6,859  x  51 
hypothetical  crystals  which  could  be  rank-ordered  by  density. 
Next,  rigid-molecule  lattice  energy  minimizations  were  per¬ 
formed  on  the  500  highest-density  candidates  using  the 
WMIN  lattice  energy  optimization  codes.26  In  this  procedure, 
the  cell  parameters,  location,  and  orientations  of  the  molecules 
are  optimized  within  space  group  symmetry  constraints  using 
the  potential  function  described  hereafter.  The  starting  values 
of  these  parameters  were  those  generated  by  the  MOLPAK 
step.  The  coordination  geometries  are  not,  in  general, 
preserved  by  WMIN  during  the  energy  minimization  and  the 
geometries  given  in  the  tables  denote  only  the  starting  configu¬ 
rations.  In  order  to  identify  duplicate  crystal  structures  pro¬ 
duced  during  the  WMIN  lattice  energy  minimization  step, 
reduced  cell  parameters  are  compared,  since  such  parameters 
provide  a  unique  representation  of  a  crystal  lattice. 

The  WMIN  lattice  energy  program  cannot,  at  this  time, 
accommodate  the  more  complicated  published  form  of  the 
SAPT(DFT)  potential.10  Rather,  WMIN  requires  that  the 
lattice  energy  be  described  by  a  sum  of  intermolecular 
atom-atom  terms  including  the  exponential,  inverse  sixth 
power,  and  Coulomb  functions: 

Vap{r)  =  Axpexp(-Bxlir)  -  Ca/s/r6  +  V^p(r),  (1) 

where 

Fc(r)=M£  (2) 
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Table  1  Parameters  for  SAPT(DFT)exp_6  potential  used  in  (he 
MOLPAK/WMIN  calculations 


Pair  (a-/l) 

Aapjkca\  mol  1 

BJ  A-1 

Cap/ A6  kcal  mol  1 

C-C 

2400957.6219 

4.7141 

331.7437 

N-C 

-380.0979 

2.0758 

266.5829 

N-N 

899731.7524 

4.7263 

214.2209 

o-c 

359683.8587 

4.3728 

322.4757 

O-N 

112623.7784 

4.2737 

259.1353 

0-0 

42389.9472 

3.7509 

313.4666 

H-C 

431554.8893 

6.0138 

172.4754 

H-N 

926.9395 

2.5410 

138.5979 

H-O 

6939.8408 

3.6332 

167.6569 

H-H 

2647.1147 

3.4678 

89.6709 

In  these  formulas,  r  is  the  interatomic  distance  between  atoms 
a  and  /?,  qa  and  qp  are  the  permanent  electric  charges  on  the 
atoms,  and  the  remaining  quantities  are  adjustable  parameters 
of  the  fit  (we  will  denote  this  potential  as  “exp-6”).  Therefore, 
we  fitted  an  exp-6  form  (denoted  hereafter  as  SAPT(DFT)exp-6) 
to  the  SAPT(DFT)  points  used  to  generate  the  SAPT(DFT) 
potential  and  used  this  function  in  crystal  structure  prediction 
methods.  The  parameters  of  SAPT(DFT)exp_6  are  given  in 
Table  1.  The  charges  are  the  same  as  in  ref.  10.  The  coefficients 
Cap  were  fitted  assuming  the  combining  rule:  Cxp  =  \/CaaCpp 
and,  therefore,  only  the  diagonal  elements  in  the  last  column 
of  Table  1  are  independent. 

NsT-MD  simulations  of  polymorphs  of  RDX 

NsT-MD  simulations  using  the  DL_POLY  version  2  suite 
of  molecular  dynamics  simulation  software27  and  the 
SAPT(DFT)  potential  of  the  RDX  dimer10  were  performed 
for  the  14  lowest-energy  crystal  structures  generated  by  the 
MOLPAK/WMIN  calculations  described  above.  The  rigid- 
molecule  approximation  was  assumed  and  symmetry  constraints 
were  not  imposed.  The  supercells  used  in  these  simulations 
were  composed  of  blocks  of  unit  cells,  with  the  contents 
initially  arranged  in  the  configuration  corresponding  to  that 
of  the  structure  produced  by  the  WMIN/SAPT(DFT)exp-6 
optimization.  The  sizes  were  selected  to  ensure  that  the 
perpendicular  widths  between  opposing  faces  of  the  simulation 
cell  were  at  least  twice  the  interaction  potential  cutoff  distance 
(12.0  A).  Coulombic  interactions  were  handled  using  Ewald 
summations.  Atomic  velocities  corresponding  to  T  =  298  K 
were  then  assigned  and  a  trajectory  of  30000  steps  (1  time 
step  =  1  fs)  was  performed  to  equilibrate  the  system.  The 
atomic  velocities  were  scaled  at  every  5th  step  during  this 
equilibration  trajectory  to  drive  the  system  to  the  desired 
temperature  more  rapidly. 

A  subsequent  NsT-MD  trajectory  consisting  of  100000 
steps  was  then  performed  in  which  thermodynamic  properties 
and  crystal  parameters  were  averaged;  atomic  configurations 
were  also  recorded  at  every  500th  step  during  the  production 
trajectory.  The  snapshots  were  used  to  generate  time-averaged 
information  about  the  locations  and  orientations  of  the  molecules 
within  the  unit  cells.  These  are  given  in  terms  of  center-of-mass 
fractional  coordinates  (sx,  sv,  sz)  and  the  Euler  angles  ( 9 ,  (j>,  t/i) 
that  transform  the  principal  axes  of  inertia  of  each  molecule  to 
its  space-fixed  coordinate  system  (i.e.,  the  Cartesian  axes).  For 
each  polymorphic  form  studied,  the  molecular  parameters  for 


each  of  the  symmetry-equivalent  molecules  in  the  unit  cell 
were  averaged  over  time  and  over  all  unit  cells  within  the 
simulation  supercell. 

Results 

Table  2  provides  the  MOLPAK/WMIN/SAPT(DFT)exp_6 
predictions  of  the  density  and  lattice  energy  (per  molecule) 
for  the  lowest-energy  crystal  for  each  of  the  51  coordination 
geometries  sampled  in  this  study  (i.e.,  only  the  global 
minimum  for  all  packing/minimizatons  starting  from  a  given 
coordination  geometry  is  given).  Many  of  the  MOLPAK 
structures  that  belong  to  the  same  space  group  but  have 
different  coordination  geometries  used  in  the  packing  portion 
of  the  calculation  (e.g.,  CB  and  CC)  converged  to  the  same 
minimum  in  the  WMIN  lattice  energy  refinement  step.  Also, 
there  are  numerous  instances  in  which  crystals  in  different 
space  groups/coordination  geometries  have  similar  lattice 
energies  and  densities;  however,  comparisons  of  reduced  cell 
parameters  for  such  cases  indicate  that  these  crystals  are  not 
duplicates.  As  evident  in  Table  2,  the  lowest-energy  crystal  (the 
same  structure  was  obtained  in  optimizations  starting  from 
coordination  geometries  CB  and  CC)  is  lower  in  energy  than 


Table  2  Density  and  lattice  energy  of  the  lowest-energy  crystal  for 
each  starting  coordination  geometry  generated  in  the  MOLPAK/ 
WMIN/SAPT(DFT)exp_6  calculations 


Space  group 

Coordination  geometry 

p/gcc  1 

U/kcal  mol  1 

Pbca 

CB.  CC 

1.860 

-31.67 

P2Jc 

AM.  FC 

1.848 

-29.73 

P2,/c 

AI 

1.814 

-29.63 

P21/C 

AK 

1.849 

-29.28 

P2\2\2\ 

AQ,  AZ 

1.751 

-28.54 

Pna2\ 

AV.  BD.  AU,  AS 

1.749 

-28.52 

C2jc 

DD 

1.748 

-27.79 

Pca2\ 

AY 

1.749 

-27.73 

P2\ 

AH,  AF 

1.758 

-27.71 

P2Jc 

FA 

1.843 

-27.68 

Pca2\ 

BH 

1.784 

-27.63 

P\ 

AB,  CA 

1.785 

-27.48 

C2/c 

DC.  DE 

1.765 

-27.28 

PI 

AA 

1.705 

-27.17 

Cc 

DA 

1.683 

-27.17 

Pbcn 

CD 

1.739 

-27.15 

P2A2 1 

AP 

1.712 

-26.98 

Pna2{ 

BF 

1.740 

-26.77 

Pbcn 

CE 

1.743 

-26.61 

P2A2 , 

BA 

1.709 

-26.53 

Pc 

AD 

1.683 

-26.31 

/>2,2,2, 

BB 

1.677 

-26.26 

Pc 

AG 

1.651 

-26.23 

C2 

DB 

1.727 

-25.89 

P2jc 

AJ,  AL.  FD 

1.711 

-25.57 

Pnn2 

AR,  AT 

1.604 

-25.36 

Pnn2 

BE 

1.631 

-25.35 

Pba2 

AW,  BG 

1.693 

-25.29 

Pnc2 

BI 

1.631 

-24.49 

Pnc2 

AX 

1.630 

-24.48 

P2\jm 

AN 

1.646 

-24.28 

P2jm 

FB 

1.646 

-24.28 

P2 

AE 

1.665 

-23.86 

P222, 

BC 

1.641 

-23.68 

Pmn2\ 

BJ 

1.508 

-23.66 

Pma2 

BK 

1.599 

-22.54 

Pm 

AC 

1.323 

-18.68 

P2m 

AO 

1.395 

-18.44 
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all  other  low-energy  crystals  for  the  remaining  coordination 
geometries  by  at  least  1.94  kcal  mol-1. 

Table  3  provides  crystallographic  parameters  (given  in  the 
conventional  setting  rather  than  reduced  cell  representations) 
for  the  14  lowest-energy  crystal  structures  generated  from 
calculations  using  the  MOLPAK/WM1N  procedure  and  the 
SAPT(DFT)exP-6  potential.  Several  structures  appearing  in 
Table  3  that  were  not  present  in  Table  2  are  subsequent  local 
minima  for  packing/minimizatons  starting  from  a  given  coordina¬ 
tion  geometry.  In  particular,  the  AM  geometry  of  the  Pl\jc  group 
contributes  a  large  number  of  low-energy  minima.  The  solutions 
are  numbered  with  increasing  energy  separately  for  each  coordina¬ 
tion  geometry  (column  3).  Also  in  Table  3,  the  time-averaged 
crystallographic  parameters  for  13  of  the  corresponding  crystals 
produced  from  NsT-MD  simulations  using  the  SAPT(DFT) 
potential  are  provided  for  comparison.  In  one  case,  the  crystal 
produced  by  the  simulation  did  not  resemble  the  starting  structure 
(MOLPAK/WMIN  solution  2,  coordination  geometry  AK,  space 
group  Pli/c),  thus  suggesting  that  this  structure  is  not  a  stable 
minimum  on  the  SAPT(DFT)  potential  energy  surface.  No  further 
analysis  regarding  this  issue  was  performed. 

In  order  to  determine  whether  space  group  symmetry  was 
preserved  during  the  NsT-MD  simulations  for  each  poly¬ 
morph,  a  unit  cell  for  a  hypothetical  ideal  RDX  crystal  was 
generated  by  applying  the  appropriate  crystalline  space  group 
operations  to  averaged  atomic  positions  for  one  of  the 
symmetry-equivalent  molecules  in  the  time-averaged  unit  cell 
produced  from  the  NsT-MD  simulation.  The  molecular  para¬ 
meters  (center-of-mass  fractional  and  orientational  Euler 
angles)  of  the  time-averaged  unit  cells  were  then  compared 
with  those  of  the  “ideal”  crystals.  Table  SI  of  the  ESIf 
provides  the  maximum  deviations  of  the  monomer’s  positions 
and  orientations  in  the  time-averaged  unit  cells  from  those  of 
the  ideal  RDX  crystals.  The  maximum  deviation  of  the  center- 
of-mass  fractional  coordinates  from  ideality  among  all  of  the 
polymorphs  is  0.0032  fractional  units  (solution  2,  AM)  and  the 
maximum  deviation  of  Euler  angles  from  ideality  is  3.01° 
(solution  1,  AK).  The  good  agreement  of  the  NsT-MD 
averaged  values  with  those  of  ideal  crystals  demonstrate  that 
the  crystalline  space  group  symmetries  are  well  maintained 
throughout  the  simulations  for  these  polymorphs. 

Deviations  of  the  NsT-MD  results  from  those  of  the 
MOLPAK/WMIN  predictions  are  also  given  in  Table  3.  The 
deviations  from  the  MOLPAK/WMIN  result  were  no  greater 
than  5.6%  in  cell  dimension  (solution  1,  AK),  1.3°  in  cell 
angle,  and  8%  in  density  (solution  1,  AK,  and  solution  1,  AV). 

A  more  rigorous  comparison  of  the  contents  of  the  cells 
generated  through  NsT-MD  simulations  and  MOLPAK/ 
WMIN  calculations  are  performed  through  superposition  of 
the  NsT-MD  time  averaged  unit  cells  onto  the  MOLPAK/ 
WMIN  counterpart  using  the  procedure  described  by  Rice  and 
Sorescu.28  The  factors  Ad  and  Ax  in  Table  3  are,  respectively, 
the  maximum  rigid-body  rotational  displacement  (in  °)  and 
translational  displacement  (in  A)  of  any  of  the  symmetry- 
equivalent  molecules  in  the  time-averaged  unit  cell  relative  to 
those  of  the  MOLPAK/WMIN  counterpart.  The  time- 
averaged  crystals  corresponding  to  solution  4,  CB,  and 
solution  1,  AK,  had  the  largest  deviations  in  rigid-body 
rotational  and  translational  displacements,  respectively. 


The  NsT-MD  lattice  energies  are  all  higher  than  the 
MOLPAK/WMIN  results  and  the  densities  are  consistently 
lower.  These  trends  can  be  partially  attributed  to  thermal 
effects  introduced  in  the  MD  simulations,  whereas  the 
MOLPAK/WMIN  calculations  represent  a  0  K  result;  thermal 
expansion  should  generate  lower  densities  and  higher  lattice 
energies.  Also,  there  are  differences  in  interaction  energies 
between  the  SAPT(DFT)exp-6  and  SAPT(DFT)  potentials:  at 
the  51  RDX  dimer  potential  energy  surface  minima  listed  in 
the  supplementary  materials  of  ref.  10,  the  root  mean  square 
deviation  from  the  ab  initio  computed  values  was  0.40  and 
0.14  kcal  mol-1,  respectively,  whereas  the  maximum  devia¬ 
tions  were  1.35  and  0.65  kcal  mol-1,  respectively.  Relative 
energy  rankings  of  the  MOLPAK/WMIN  and  NsT-MD 
structures  also  differ,  as  is  apparent  in  Fig.  1 . 

Lattice  parameters  and  density  for  the  experimental  crystal 
are  given  in  Table  3  for  comparison  with  both  the  Pbca  low- 
energy  structure  from  the  MOLPAK/WMIN  calculations  and 
the  corresponding  time-averaged  structure  generated  in  the 
NsT-MD  simulations.  Agreement  with  experiment  for  all  cell 
parameters  is  very  good,  especially  those  of  the  NsT-MD 
results.  The  deviations  from  experiment  of  the  NsT-MD 
densities  and  cell  edge  lengths  in  the  a,  b  and  c  directions 
are  —0.8,  0.4,  0.3  and  0.1%,  respectively.  Deviations 
from  experiment  in  the  MOLPAK/WMIN  case  are  slightly 
larger  (3.0%,  —0.7,  —0.8  and  —1.5%,  respectively),  partially 
due  to  lack  of  inclusion  of  thermal  effects.  Table  S2  of  the 
ESI|  provides  a  comparison  of  time-averaged  center-of-mass 
fractional  coordinates  and  orientational  Euler  angles  with 
the  corresponding  experimental  values  for  RDX  at  298  K 
and  1  atm.  Agreement  with  experiment  for  these  para¬ 
meters  is  also  outstanding,  with  the  largest  deviation  in 
molecular  orientation  occurring  for  the  Euler  angle  i/r,  which 
differs  from  experiment  by  ~2.7°.  Also,  the  largest  deviation 
of  the  location  of  the  molecular  mass  center  is  less  than  0.07  A. 
This  excellent  agreement  is  visualized  in  Fig.  1  of  ref.  9. 

Further  tests  of  the  SAPT(DFT)  potential  were  performed 
to  assess  its  ability  to  predict  the  temperature  and  pressure 
dependences  of  the  RDX  crystal.  The  simulation  cells  used  in 
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Fig.  1  Lattice  energies  (kcal  mol  ')  for  the  14  lowest-energy  confor- 
mers  generated  by  MOLPAK/WMIN.  Solid  circles  denote  the  results 
from  MOLPAK/WMIN  (0  K)  calculations  using  SAPT(DFT)exp_6  and 
the  open  circles  denote  the  results  from  NsT-MD  (298  K,  1  atm) 
simulations  using  the  SAPT(DFT)  potential. 
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Table  4 

Temperature-dependent  lattice  parameters  and  thermal  expansion  coefficients  (y)° 

Tj  K 

a/A 

b/A 

c/A 

a/° 

f!/° 

y/° 

v/A3 

100 

13.1829 

11.5526 

10.6758 

90.00 

90.00 

90.00 

1625.893 

200 

13.2207 

11.5915 

10.7133 

90.00 

90.00 

90.00 

1641.783 

250 

13.2385 

11.6139 

10.7333 

90.00 

90.00 

90.00 

1650.244 

273.15 

13.2487 

11.6225 

10.7439 

90.00 

90.00 

90.00 

1654.368 

298 

13.2591 

11.6338 

10.7542 

90.01 

90.00 

89.99 

1658.873 

325 

13.2693 

11.6459 

10.7670 

90.00 

90.01 

90.00 

1663.857 

28.00  (26.41)  [34.8] 

33.91  (86.96)  [47.5] 

34.59  (79.26)  [47.2] 

95.42  (192.6)  [129.6] 

a  Units  for  linear  and  volume  expansion  coefficients  are  in  10 

6  K  .  Values  in  parentheses  are  derived  from  experimental  measurements  given  in 

ref.  29,  and  values  in  brackets  are 

theoretical  values  of  Sorescu  et  al.30 

Table  5 

Lattice  parameters 

as  functions  of  ] 

pressure 

P/ GPa 

a/A 

b/A 

c/A 

v/A3 

«/°  Pf 

NsT-MDc 

y/° 

Expt" 

NsT-MD* 

Expt" 

NsT-MD* 

ExpU 

NsT-MD* 

Expt" 

NsT-MD* 

0.00 

13.20 

13.2591 

(0.4) 

11.60 

11.6338 

(0.3) 

10.72 

10.7542 

(0.3) 

1641.446 

1658.873 

(1.1) 

90.01 

90.00 

89.99 

0.73 

13.02 

13.0701 

(0.4) 

11.42 

11.4568 

(0.3) 

10.53 

10.5561 

(0.2) 

1564.298 

1580.697 

(1.0) 

90.00 

90.00 

90.00 

1.75 

12.89 

12.9089 

(0.1) 

11.22 

11.2925 

(0.6) 

10.27 

10.3735 

(LO) 

1485.509 

1512.182 

(1.8) 

90.00 

90.00 

90.00 

2.75 

12.80 

12.8004 

(0.0) 

11.07 

11.1721 

(0.9) 

10.16 

10.2400 

(0.8) 

1441.190 

1464.399 

(1.6) 

90.00 

90.00 

90.00 

3.36 

12.74 

12.7484 

(0.1) 

11.00 

11.1106 

(1.0) 

10.08 

10.1701 

(0.9) 

1413.285 

1440.509 

(1.9) 

90.00 

90.00 

90.00 

3.95 

12.67 

12.7055 

(0.3) 

10.93 

11.0568 

(1.2) 

10.03 

10.1083 

(0.8) 

1388.664 

1420.041 

(2.3) 

90.00 

90.00 

90.00 

"  Ref.  31 

.  *  Values 

>  in  parentheses 

are  percent 

deviations  from  experimental  values. 

£'  Only  predicted  cell  angles  are 

reported;  experimental  values 

are  all  90°. 


the  calculations  at  various  temperatures  and  pressures  were 
composed  of  3  x  3  x  3  crystallographic  unit  cells  and  were 
allowed  to  equilibrate  to  the  specified  thermodynamic  condi¬ 
tion  before  properties  were  calculated  and  accumulated  for 
averaging.  First,  a  series  of  NsT-MD  calculations  were  per¬ 
formed  over  the  temperature  range  100-325  K.  Temperature- 
dependent  crystallographic  parameters  are  given  in  Table  4. 
Linear  and  volume  coefficients  of  thermal  expansion  (CTE) 
are  extracted  from  these  data  and  are  also  given  in  Table  4 
along  with  experimental29  and  other  theoretical  values30  for 
comparison.  We  note  that  the  Sorescu  et  al.  study30  reported 
experimental  linear  coefficients  that  are  one-third  of  the 
reported  measured  volume  CTE  value;  however,  the  experi¬ 
mental  information  indicates  that  thermal  expansion  of  the 
lattice  is  highly  anisotropic.  Both  the  work  of  Sorescu  et  al. 
and  our  predictions  here  indicate  anisotropy  in  the  thermal 
expansion,  but  to  a  significantly  lesser  extent  to  that  observed 
experimentally.  The  predicted  CTE  for  the  a  axis  is  in  excellent 
agreement  with  the  measured  value;  however,  the  CTEs  for  the 
remaining  two  axes  (and  subsequently  the  volume  CTE)  are 
too  low  by  factors  of  2-2.5.  It  is  possible  that  some  of  the 
discrepancies  could  be  due  to  the  rigid-molecule  approxima¬ 
tion  imposed  in  the  NsT-MD  simulations. 

The  pressure  dependences  of  the  lattice  parameters  were 
also  calculated  for  pressures  ranging  from  0  to  3.95  GPa,  and 
are  given  in  Table  5  along  with  experimental  values.  As  shown 
in  the  table,  predicted  values  in  the  edge  lengths  deviate  from 
experiment  by  no  more  than  1.2%,  and  the  deviation  of  the 
volumes  over  the  pressure  range  varies  from  1-2.3%.  A  visual 
comparison  with  experimental  values31  is  given  in  Fig.  2.  The 
predicted  data  were  fitted  to  the  Murnaghan  equation  of  state 


Pressure  (GPa) 

Fig.  2  Variation  of  predicted  and  experimental  unit  cell  volume  and 
lattice  dimensions  with  pressure  for  RDX.  Experimental  values  are 
taken  from  ref.  3 1 . 
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to  provide  estimates  of  the  bulk  modulus  B0  at  zero  pressure 
and  its  pressure  derivative  BQ'.  The  predicted  values  of  B0  and 
B0'  (12.7  GPa  and  8.2,  respectively)  are  in  good  agreement 
with  reported  values  derived  from  experimental  measurements 
(13.0  GPa  and  6.6,  respectively).31 

4.  Summary 

Crystal  structure  prediction  methods  employing  the  MOLPAK/ 
WMIN  suite  of  codes  were  used  to  explore  various  poly¬ 
morphs  of  the  RDX  crystal  on  the  SAPT(DFT)exp_6  potential 
energy  surface.  The  low-energy  crystal  produced  by  these 
calculations  matched  the  experimental  crystal  at  ambient 
conditions.  This  structure,  as  well  as  the  13  other  lowest- 
energy  crystals  generated  by  these  calculations,  were  used  as 
initial  structures  in  NsT-MD  simulations  at  T  =  298  K  and 
1  atm  pressure.  The  NsT-MD  simulations  used  the  more 
elaborate  SAPT(DFT)  potential  energy  description.10  For  all 
but  one  of  the  structures,  the  NsT-MD  predictions  of  crystallo¬ 
graphic  parameters  are  in  good  agreement  with  the 
MOLPAK/WMIN  predictions,  with  a  maximum  5.6%  devia¬ 
tion  in  cell  lengths,  1.3°  in  cell  angles  and  8%  in  densities. 
Additionally,  the  space  group  symmetries  are  maintained 
throughout  the  NsT-MD  simulations,  except  for  one  case. 
The  NsT-MD  results  show  a  slight  decrease  in  density  and 
higher  lattice  energies,  which  can  be  attributed  to  thermal 
effects  and  differences  in  the  potential  energy  functions.  Also, 
the  relative  energy  rankings  differ  between  the  two  models; 
however,  the  low-energy  structures  for  both  sets  of  calcula¬ 
tions  correspond  to  the  experimental  crystal.  The  crystallo¬ 
graphic  parameters  and  contents  of  the  low-energy  crystal 
produced  from  the  NsT-MD  simulations  at  room  conditions 
matched  those  of  the  experimental  crystal  to  within  1  %  of  cell 
edge  lengths  and  density  and  0.01°  of  cell  angles.  Further  tests 
of  the  SAPT(DFT)  potential  explored  its  ability  to  predict  the 
temperature  and  pressure  dependences  of  the  crystallographic 
parameters  of  RDX.  The  results  indicated  that  changes  in 
crystallographic  parameters  with  pressure  are  similar  to 
experimental  values;  additionally,  the  bulk  modulus  and  its 
pressure  derivative  are  in  good  agreement  with  experimentally 
derived  values.  Predicted  thermal  expansion  coefficients  indi¬ 
cated  that  expansion  of  the  lattice  is  only  slightly  anisotropic, 
in  contrast  with  experiment.  Also,  only  one  of  the  linear 
expansion  coefficients  was  in  good  agreement  with  experi¬ 
mental  values;  the  remaining  two  coefficients  were  sub¬ 
stantially  smaller  than  the  experimentally  derived  values. 
However,  it  is  possible  that  the  thermal  expansion  behavior 
might  be  strongly  influenced  by  intramolecular  motion,  which 
was  not  present  in  our  simulations.  The  ability  of  the  model  to 
predict  the  pressure  dependence  of  the  lattice  parameters  is 
much  better,  with  a  maximum  deviation  of  2.3%  in  the  volume 
and  1 .2%  in  cell  lengths  over  the  entire  pressure  range  that  was 
explored.  The  success  of  the  SAPT(DFT)  potential  in  describ¬ 
ing  condensed  phase  RDX  at  the  ambient  state  and  com¬ 
pressed  states  provides  great  incentive  to  use  this  approach  in 
developing  such  analytical  potential  energy  functions  for  use 
in  atomistic  simulations  of  the  condensed  phase  of  other 
molecular  crystals. 
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